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Elastic systems that are spatially heterogeneous in their mechanical response pose special challenges for 
molecular simulations. Standard methods for sampling thermal fluctuations of a system's size and shape proceed 
through a series of homogeneous deformations, whose magnitudes can be severely restricted by its stiffest parts. 
Here we present a Monte Carlo algorithm designed to circumvent this difficulty, which can be prohibitive in 
many systems of modern interest. By deforming randomly selected subvolumes alone, it naturally distributes 
the amplitude of spontaneous elastic fluctuations according to intrinsic heterogeneity. We describe in detail 
implementations of such "slice moves" that are consistent with detailed balance. Their practical application is 
illustrated for crystals of 2D hard disks and random networks of cross-linked polymers. 



I. INTRODUCTION 

Modern intersections of chemistry, biology, and materi- 
als science focus attention on systems that are substantially 
nonuniform in their spatial organization; examples include in- 
terfaces, structural elements of the cell such as the cytoskele- 
ton, and systems undergoing phase transitions. Tools of sta- 
tistical mechanics that could help clarify their structure and 
function often do not apply straightforwardly or efficiently in 
the face of such heterogeneity. This paper concerns a class of 
computational methods that suffer in this way. 

Specifically, we address methods for simulating shape fluc- 
tuations of elastic materials. Pioneered by Parrinello and 
Rahmanf]]] in the context of molecular dynamics, these ap- 
proaches extended constant pressure simulation techniques 
andlll [H, like their predecessors, opened the door for novel 
computational studies of phase transitions^ H, S 0]. The 
basic idea of these approaches is simple to understand: treat 
the parameters of a system's overall geometry as fluctuating 
dynamical variables, on the same footing as molecular coor- 
dinates. In practice, it is convenient to isolate changes in box 
size and shape by introducing scaled (reduced) coordinates, 
h = ^h^ r i (using Einstein summation convention), where r'- 
is the j-coordinate of the position vector of atom i and hij is 
a matrix of d lattice vectors defining the periodically repli- 
cated (i-dimensional box geometry. Parrinello and Rahman 
constructed a Lagrangian with fictitious terms involving h(j 
and its time derivatives, allowing dynamical simulations of a 
system with fluctuating shape. One can similarly use Monte 
Carlo simulations to sample the components of hu from a 
Boltzmann distribution[2|]. 

The problem with applying these methods to heterogeneous 
materials is also simple to understand. When hij changes, so 
do the physical positions of all atomic coordinates. For ex- 
ample, if one of the basis vectors in defining a rectangular 
simulation cell is scaled by some factor, the corresponding 
components of all position vectors r, become multiplied by 
the same factor. The stiffness of the resulting motion is de- 
termined by the resistance of molecular interactions to these 
scale and shear transformations. Sampling efficiency is thus 
determined by the proverbial "weakest link": If even a small 
part of the system strongly resists deformation, then a sim- 
ulation must await rare, transiently softening fluctuations in 



local structure that facilitate changes in overall geometry. A 
system featuring many locally stiff regions becomes nearly in- 
tractable, since the likelihood of many rare local fluctuations 
occurring simultaneously is extremely small. 

Our attention to these methodological issues is driven by 
an interest in the polymer networks that determine elastic 
properties of living cells |8t |9|]. As a crude but illustrative 
model of these cytoskeletal materials, consider a collection 
of semi-flexible filaments, placed and oriented at random on 
a two-dimensional plane, that are permanently cross-linked 
wherever they intersect OlOll . In this case spatial variations in 
cross-link density effect substantial variations in local stiff- 
ness. Even modest global strains are typically not tolerated by 
the densest regions of the network. A Monte Carlo simulation 
of such a model could achieve a reasonable acceptance rate 
only by using very small displacements in system geometry. 
As a result, relaxation would proceed quite sluggishly. 

In this paper we present a technique that can remove 
these difficulties by allowing heterogeneous deformations. By 
transforming only part of a system, we avoid hinging fluc- 
tuations of the system as a whole on its stiffest parts. We 
consider such motions as trial moves in a Metropolis Monte 
Carlo scheme. We term these as "slice moves", since they 
proceed by choosing slices of a system that deform, leaving 
the remainder of the system internally unaffected. In section 
HI1 we introduce the method in detail for both constant pres- 
sure (NPT) simulations and constant stress simulations, pay- 
ing careful attention to the requirement of detailed balance. 
We illustrate the method in section ITlTl through application to 
the elasticity of crystals of 2-dimensional hard disks, and ran- 
dom networks of cross-linked semi-flexible polymers, and in 
section [TVl we conclude. 



II. FAST SAMPLING THROUGH PARTIAL VOLUME 
MOVES 

The basic flaw of conventional strain sampling techniques, 
when applied to nonuniform systems, is their global na- 
ture. We localize strain moves in Monte Carlo simulations 
by choosing thin slices of a system, outside of which inter- 
molecular geometries are undisturbed. Fig. Q] illustrates such 
a partial volume move. In this two-dimensional example, 
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FIG. 1: A partial volume move. The area v transforms to v', co- 
transforming the shaded areas. The rest of the simulation box re- 
mains unchanged; the new periodically replicating simulation box 
boundaries are shown as dotted lines. 



subvolumes to be deformed are defined by two intersecting 
swaths. As a trial move, we deform the region v shared by 
both slices, producing a new subvolume geometry v'. The 
requirement that regions outside the two slices remain unde- 
formed then uniquely determines transformations within the 
remaining slice regions (i.e., within one but not both swaths). 
By choosing the slices' locations and widths at random, we 
can in effect sample around problematically rigid parts of a 
configuration. 

Algorithmically, such a "slice move" proceeds as follows: 

1 . Select a particle at random, whose position rW serves 
as an anchor for the primary deformation subvolume v. 

2. Select random parallelepipeds v and v' defining initial 
and final geometries of the subvolume. 

3. Determine additional parallelepipeds v^ l \ V^ 2 \ etc., that 
connect v with its periodic images (See Fig. |2j. These 
regions, together with v, form the intersecting slices that 
will be deformed. Repeat for v'. 

4. Calculate each particle's position in the deformed trial 
state according to strains applied to the region in which 
it resides. 

5. Evaluate the change in internal energy AU (accounting 
for the change in periodic boundary conditions) and the 
work W associated with external forces. 

6. Accept or reject the trial move with a probability deter- 
mined by the total change in energy relative to T 

We will describe two variants of such a trial move. The 
simpler version involves only the limited class of transforma- 
tions that switch between rectangular system geometries, for 
which the shape matrices describing v, v', h, etc. are all diag- 
onal. The more general, and in practice much more compli- 
cated version, includes the possibility of shear deformations 
as well. 

Different periodic images of a particle may move differ- 
ently in the course of a partial volume move. Detailing the 
algorithm is therefore greatly simplified by a careful and spe- 
cific choice of images. Fig.|2]illustrates how we select among 
each particle's set of periodically replicated coordinates, ac- 
cording to the subvolume it occupies. 
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FIG. 2: The connecting subvolumes v' 1 ) and v^ 2 ), and the coordinate 
shifts and periodic boundary conditions during a slice move. The 
dashed lines show the original periodic boundary box outlined by hjj 
(already shifted to r' ' at its origin), while the base periodic image 
during the slice move is shown with the thick lines. 



First, we require that each subvolume (v, v' 1 ', v^ 2 \ and the 
unperturbed region u) is not fragmented across system bound- 
aries. Since the subvolumes are themselves repeated in space, 
this criterion does not by itself uniquely specify a choice of 
particle images. We further choose that the un-fragmented 
regions are adjacent in a particular way: the subvolumes v^ 1 ), 
v^ 2 ), and u must all contact the anchor point A ' . Note that this 
scheme results in a collection of particle coordinates that do 
not lie within the boundary of a single simulation cell. Finally, 
we translate all particles uniformly so that the anchor point 
A " 1 lies at the origin. While geometrically straightforward, 
this set of operations carries a nontrivial computational over- 
head. We describe an efficient implementation in AppendixlAl 



A. Scaling slice moves in a rectangular simulation box 

The simplest slice move modifies only the scale of rectan- 
gular slices along the corresponding lattice vectors, as shown 
in Fig. [3] Such a move effects a change in system volume and 
aspect ratio, but does not change the relative directions of lat- 
tice vectors defining periodic boundary conditions. Here, we 
will take the box matrix /i,-, to be purely diagonal, both before 
and after the distortion. B21I1 

The width of initial and final slices, together with the slice 
origin, completely specify an instance of this partial volume 
move. We define v, as the length of subvolume v in direction i, 
and vj as its length in the trial configuration. The deformation 
Si is simply determined by the ratio of these widths, 



(1) 



Recall that the slice origin r^°- is assigned to be the location 
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FIG. 3: A partial scaling move in a rectangular box. For notational 
simplicity we label co-transforming regions Vj along directions j as 
v x (for j = 1) an d v v (f° r j = 2). Corresponding regions in the de- 
formed system are similarly labeled v' x and vi rather than Vj and v' 2 , 
respectively. 



of a randomly selected particle. By construction, the lengths 
Ui of the undeformed region u do not change during a slice 
move. The box matrix for the trial configuration is therefore 
given by 



h'ij = (at + v'i) dij 



(2) 



Accounting for this deformation, the reference frame trans- 
lation placing A°> at the origin, and the choice of periodic im- 
ages depicted in Fig.|2j we can write the position of particle i 
in the trial configuration as 



s^i if n > o, 

r, otherwise. 



(3) 



In practice, the partial volume transformation is more con- 
veniently performed using reduced coordinates . In this rep- 
resentation the coordinates of particles in the unperturbed re- 
gion u change even though their physical arrangements do not: 



Ahy rj otherwise. 



(4) 



Here, the relative deformation matrices A/z|' nt ' and A/z,-; are 



given by 



Ah 



(int) 



: h' ik l s k h kj , 



(5a) 
(5b) 



In a Metropolis Monte Carlo simulation, the probability 
^acc(r — ► r') with which a partial volume move from mi- 
crostate F to microstate T' should be accepted is dictated by 
the requirement of detailed balance: 



Pacc(r^r') 



PequnPgenOT') 



' Pequ(r) Pgen(r'|r) 



(6) 



Here, p eq u(r) is the equilibrium weight of microstate F in the 
thermal ensemble of interest; and P 2en (r'|r) is the conditional 
probability distribution for generated trial configurations F', 



given the original configuration F. This generation probabil- 
ity depends on the way in which slice geometries are chosen. 

Let 77 (v mi v' nl r^^j be the distribution of parameters specify- 
ing a partial volume move. Since the resulting microstate is 
uniquely defined by Eqs.[Tll5bl P ge n(r'|r) can be written as a 

product of 7} n>, v',A ^ and Dirac delta functions describing 
the coordinate transformations 



p g en(r'|r) = nn 5 W 

P=l {i)„ 
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(7) 



The notation JT_(i) indicates a product over all N p particles 
whose positions are influenced by partial volume scaling in 
the direction p, i.e., particles lying within the slice that runs 
perpendicular to the p th lattice vector. Similarly, Yl[j] p de- 
notes a product over the — N p particles whose coordinates 
are unaffected by scaling in the p direction. Note that the an- 
chor point does not change during the transformation - for 
accounting purposes, the corresponding particle lies outside 
the deformed sub volumes. 

If the distribution rj (v, v',r' ^ is symmetric with respect 

to exchange of v and v' (i.e., if the original subvolume v and 
the distorted subvolume v' are selected in the same way), the 
ratio of generation probabilities appearing in Eq. [6] evaluates 
simply to 



/Wry) 
^gen(r'ir) 



[l s r 



(8) 



For a system held at fixed temperature T = (ksfi)^ 1 and 
isotropic pressure P, equilibrium probabilities depend on in- 
ternal energy E as well as the total volume V, p e qu(r) °= 
exp[— P(E +PV)]. The corresponding Metropolis acceptance 
probability for a partial volume move is then 



Paec(r^r') = 



min 




:exp(-j3[A£ + p(|^ r |-|v|)]) 



(9) 



where AE is the change in internal energy resulting from the 
trial deformation. |A, ; | denotes the determinant of a matrix 
A, so that \hij\ and represent the volumes of original and 
trial states, respectively. 

A partial volume move closely resembles a conventional 
global strain move when the deformation subvolume v encom- 
passes the whole system, v; = /i,;. In this case all particle co- 
ordinates (except those of the slice anchor point) are subjected 
to scaling in each direction, N p = N — 1 . The acceptance prob- 
ability then becomes min[l,e~^ A£/ ], with an effective poten- 
tial U = E + PV + (N - l)k s TlnV, much as in a standard 
isothermal-isobaric Monte Carlo simulation ! 1 ill . 
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FIG. 4: A partial volume move with shear components. The thick 
lines represent the box boundary as used in the slice move. 



B. Slice moves with shear 

Because lattice vector orientations are invariant under the 
deformations described in the preceding section, those trial 
moves do not suffice for simulating shear fluctuations. In this 
section we present a generalization of slice moves suitable for 
that purpose. It is tempting to proceed by selecting rectan- 
gular slices, as before, and then distorting them into paral- 
lelotope shapes (like the deformation sketched in Fig. [TJ. If 
restricted to rectangular slices of the initial state, however, 
a move of this sort is irreversible and therefore inconsistent 
with detailed balance. Incorporating shear correctly requires 
the possibility that slices of the initial state also be shaped as 
parallelotopes (parallelograms in d — 2 and parallelepipeds in 
d = 3), as shown in Fig. [4] 

The product of a slice move including shear components 
is a box matrix whose component vectors differ in direction 
from those of the initial state. We can therefore no longer 
treat distortions in different directions as independent defor- 
mations. As a mathematical consequence, we require matri- 
ces (rather than vectors as in the previous section) to describe 
subvolume shapes . 

Let vy be ad x d matrix whose rows are vectors spanning 
the edges of the deformation subvolume v. Similarly, the rows 
of Uij span the edges of the undisturbed region u. As sketched 
in Fig. Aperiodic boundary conditions demand that 



(10) 



Particles residing neither in v nor in u belong to one of sev- 
eral co-transforming subvolumes, whose shape matrices com- 
bine one or more rows of vy with one or more rows of My. 
For the case d — 2 we denote the two co-transforming regions 
v^ 1 ) and v^ 2 \ as shown in Fig. |4] (We will discuss the three- 
dimensional case later.) Subvolume vW connects the right 
edge of v with the left edge of its horizontally replicated peri- 
odic image; v' 2 ) connects top and bottom edges of vertically 

(k) 

replicated periodic images. Matrices Vy describing these re- 
gions share one row with vy and one row with My, 



..(*) 



Vy if/** 
My otherwise. 



(11) 



and 



h'fj = Uij + Vy 

= h U + ( v 'ij - v ij) 
/(*) _ / v'ij 



'' 1 My otherwise 



(12) 



(13) 



With these definitions we can compactly express deforma- 
tion matrices describing the strain applied to each subvolume: 



< s ij = v' ik v k j, 



(14a) 

4?=vf(v W )*/. d4b) 



Particle positions in the trial microstate can finally be written 

(15) 



Sijrj if r lies in v 

(*) 
i j 



r'j = ^ s\j rj if r lies in vW 
r, if r lies in u 



The generation probability for slice moves including shear 
is similar to that of the simpler deformations described by 
Eq.H 



p ge n(r'|r) = 77(v,vV(°>)n 
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(16) 



where S?{c£) denotes the set of N a particles that reside in 
subvolume a. If the subvolumes v and v' are selected inde- 
pendently from the same distribution, as we assumed in Eq.[8] 
then the ratio of backward and forward probabilities becomes: 



^gen(r'|r) 





(2) 


N v (2) | | 
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v%(v { \j ' vg(vW) 



\{2) 



V 'ik V k} 



N v 

(17) 



Detailed balance can therefore be satisfied by accepting these 
slice moves with a probability: 



^acc(vy,v^,ri 0) ) =min 



(t) N v m (2) ^(2) | , Nv 



x exp 



We employ similar definitions for the trial configuration, so 



AE + W„ 



(18) 
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As in Eq. [9] AE denotes the change in internal energy result- 
ing from the trial move. The mechanical work W ext against 
external forces may depend on the box matrices hjj and h'^ in 
a complicated way if applied stresses are anisotropic. For the 
simplest case of constant applied isotropic pressure, this en- 
ergy takes the familiar form of pressure-volume work, W ex t = 

ni%l-N)-@] 

Slice moves in three dimensions require a larger and 
slightly more complicated set of co-transforming subvolumes. 
We denote these six regions v^ kl \ where k and / take on inte- 
ger values corresponding to the three cardinal directions, and 
v i lk ) refers to the same region as v^ l > . The parallelepiped v'**) 
connects a face of the primary deformation subvolume v with 
the opposing face of its periodic image in direction k, much 
as for the d = 2 case. In d = 3 these subvolumes must them- 
selves be connected by co-transforming regions in order to 
preserve the undisturbed parallelepiped u. The region v^, 
for example, connects periodic images of V" in the direction 
/ (or, equivalently, periodic images of in the direction k). 
Shape matrices for these subvolumes are given by 



employ a wide range of subvolume shapes in order to accom- 
modate elastic inhomogeneities that are a priori unknown; at 
the same time, typical acceptance probabilities can be very 
low if v and v' differ substantially. Below we describe a pro- 
cedure for choosing deformation regions that addresses both 
of these goals, while respecting the necessity of statistical in- 
dependence. 

A natural method for generating random shape matrices 
would draw elements from a uniform distribution limited to 
a certain range e. This approach pits the above goals against 
one another. Small values of e discourage generating diverse 
subvolume shapes. Large values of e permit significant dis- 
parity between independent samples. One simple way of cir- 
cumventing this dilemma is to vary at random the mean values 
of distributions from which matrix elements are selected. 

Toward this end we define a symmetric reference matrix 



rnd(v? 



V,2 



rnd(- 
rnd(vi 



~mm ~max 

XV ' XV 



mm 

yy ' 



(23) 



(ki) 
v- ■ 



Vij if j + k and j + I 
Uj ; otherwise, 



(19) 



Aside from this enlarged set of subvolumes, coordinate 
transformations and generation probabilities proceed just as 
for d — 2. For example, the deformation matrix for region 



is given by 



-(H) 



im \ J 



-1 
mj 



(20) 



whose elements change stochastically over the course of a 
Monte Carlo simulation. Here, md(a,b) denotes a random 
number uniformly distributed between a and b. We employ a 
given realization of v,- ; as a random offset for selecting both 



vij and v'if. 



-Avy = ' 



'mdi-Exx,^) md(-e„,e yx ) 
Avi2 rnd(-%,%) 



(24a) 



The acceptance probability dictated by detailed balance is also 
simply generalized: 



^acc(r^r , )=min 



n 



(a) 



x exp 



AE + Wv 



(21) 



where the product runs over all subvolumes a (including v, 
m, and the co-transforming regions yW) with corresponding 
deformation matrices s,f . 



C. Selecting v and v' 

We have shown that detailed balance is straightforward to 
achieve with slice moves, provided the selection of subvol- 
umes v and v' is symmetric: 



J7 (v,v',r( )) = T7(v',v,/°)). 



(22) 



/ md(-e XXl e xx ) rnd(-e XJ ,,e^) 



Av', 
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rnd(-e v> ,,%) 



(24b) 



By controlling the ranges e^, e„ and e yy of variations about 
the reference geometry, similarity of v and v' can be assured 
and a reasonable acceptance probability maintained. Note that 
matrix symmetry allows only three elements of v i; - to be cho- 
sen independently. 

It can be demonstrated that this scheme obeys detailed bal- 
ance for any set of fixed parameters v™ n , v™ ax , and £y, so long 
as slices do not exceed the overall system size (just as con- 
ventional constant pressure simulations require volume incre- 
ments smaller than the system's total volume). This constraint 
should not be limiting: if elastic heterogeneity calls for slice 
moves, they will be useful only if typical slices are smaller 
than natural correlation lengths for strain fluctuations. 



III. SIMULATIONS 



Eq.|22]is most easily satisfied by choosing the corresponding 
shape matrices independently, and from the same distribution. 
Consequently, one's choice of deformed geometry vj ; cannot 
be biased by the system's current shape hij. This restriction 
poses a challenge to efficient sampling. It is advantageous to 



We have implemented slice moves in computer simulations 
of two model systems, both to verify that equilibrium ensem- 
bles of spontaneous box deformations are correctly sampled 
and to demonstrate improved efficiency for elastically hetero- 
geneous systems. 
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FIG. 5: Probability distribution function P(h xv ) for the shear compo- 
nent h xy of the box shape matrix for the two-dimensional hard disk 
system at pressure P = 10 k# T /a 2 . 



A. Validation: Hard disk solids in two dimensions 

The elastic properties of two-dimensional crystals compris- 
ing hard disks have been calculated with high precision!!, 
\H H Elin efforts to assess the possibility of a KTHNY 
transition |[l5i MM . (A sufficiently low Young's modulus sig- 
nals instability to the creation of dislocations, implying a two- 
stage, continuous transition to the disordered fluid phase.) 
Here we use those results, obtained using conventional ap- 
proaches, as benchmarks for validating our new methods. 

We have simulated systems of N = 780 hard disks at con- 
stant pressure P with periodic boundary conditions. Pairwise 
interactions forbid interparticle separations smaller than the 
particle diameter a but otherwise do not bias spatial arrange- 
ments. For P^t9 ksT/a 2 the equilibrium state is a crys- 
talline solid[17], whose bulk modulus B and effective shear 
constant /x e ff we determine from distributions of spontaneous 
fluctuations in Lagrangian strain^ S 0, EH]. Fig. [5] illus- 
trates the breadth of strain fluctuations for the specific case 
P = 10 k B T /a 2 . 

We have calculated B and jii e ff for hard disk solids at several 
pressures. In each case we performed one simulation using 
slice moves and one using exclusively conventional methods. 
Our results are shown in Table IHI Al along with previously 
reported values for systems at similar conditions. Simulations 
with and without slice moves agree well, yielding results for 
most pressures that lie well within error margins. At P = 10 
and P = 1 1 there are some differences in the bulk modulus B 
that might be explained through the moderate softness of the 
crystal itself; we speculate that the slice moves enable better 
sampling of local defects that lead to elastic heterogeneities, 
leading to slightly different bulk moduli II 1411 . and unrealisti- 
cally low error estimates in the case of no slice moves. Aside 
from that, the results match the most accurate results pub- 
lished elsewhere. 

In terms of computational efficiency, hard-disk crystals rep- 
resent something of a worst-case scenario for slice moves: the 
interparticle potential requires little numerical effort to evalu- 
ate, elastic response is spatially uniform, and the interactions 
are extremely short-ranged. These factors render significant 
the added computational overhead of a slice move relative to a 



method 


N 


P 


B 


Meff 


sf, slice moves 


780 


9.545 


42(1) 


15.3(3) 


sf, no slice moves 


780 


9.545 


42.4(5) 


15.0(2) 


sf, slice moves 


780 


10.0 


48.4(7) 


21.96(9) 


sf, no slice moves 


780 


10.0 


51.0(2) 


21.78(5) 


sf [12 


896 


10.0 


49.2(8) 


21.9(3) 


sf, slice moves 


780 


11.0 


60.8(8) 


26.3(1) 


sf, no slice moves 


780 


11.0 


57.9(3) 


26.1(2) 


sf, slice moves 


780 


11.6 


69(1) 


28.89(9) 


sf, no slice moves 


780 


11.6 


67(1) 


28.9(3) 


sf [12] 


896 


11.6 


68(2) 


28.8(4) 


stress-strain r 121 


7020 


11.6 


67(1) 


28.8(3) 


sf, slice moves 


780 


13.0 


85(1) 


35.8(2) 


sf, no slice moves 


780 


13.0 


85(1) 


35.8(3) 


sf, slice moves 


780 


15.4 


118(2) 


48.5(3) 


sf, no slice moves 


780 


15.4 


119(2) 


49.4(3) 


sf [12] 


896 


15.4 


118(1) 


48.7(3) 


stress-strain 1121 


7020 


15.4 


118(1) 


49.2(3) 


sf, slice moves 


780 23.1 


255(4) 


104.9(5) 


sf, no slice moves 


780 23.1 


260(5) 


104.0(8) 


sf [12 


896 


23.1 


251(4) 


104(1) 


stress-strain H21 


7020 23.1 


252(3) 


103.3(9) 



TABLE I: Elastic constants of 2D hard disk crystals, comparing 
strain fluctuation ('sf') simulations with and without slice moves, 
and previous work. N is the number of particles in the system, P 
is the applied isotropic pressure, B is the bulk modulus, and /i e ff is 
the effective shear modulus! 13]. Units are in particle diameter o" and 
k^T; the number in parentheses represents the standard error in the 
last digit. 



conventional strain move. Specifically, the run time is approx- 
imately 25% longer for simulations with slice moves com- 
pared to those lacking slice moves but comprising the same 
total number of strain moves. 



B. Cytoskeletal networks 

Slice moves may offer considerable computational savings 
whenever the resistance of a system to strain varies signifi- 
cantly in space. Here we demonstrate their utility for a model 
elastic gel inspired by the polymeric framework of living cells. 
This two-dimensional system comprises a collection of semi- 
flexible filaments connected by cross-links. For the specific 
model we consider here, cross-links enforce overlap of two 
filaments at fixed points along their contours but do not con- 
strain the angle at which they intersect. A thorough exami- 
nation of this model's elastic response will be presented in a 
forthcoming paper. 

We construct a particular realization of the network by lay- 
ing down straight filaments of fixed length, located and ori- 
ented at random, until a desired density is achieved. Wherever 
filaments intersect, they become permanently cross-linked. 
We will focus on two such configurations, differing in den- 
sity. Both are shown in Fig. [6] The contour length of a fil- 
ament segment between two-cross-links is set such that the 
initial distance between cross-links minimizes the segment's 
free energy. 
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FIG. 6: Two randomly laid down networks of model actin filaments 
during a Monte Carlo simulation. The line widths denote instanta- 
neous parallel strain for the segment, while the darkest shades de- 
note higher instantaneous perpendicular strain for the segment. The 
sizes of the systems are 0.25//, x 0.25//, ~ 2.5 flm x 2.5 flm, where 
l p is the persistence length of actin. The system in (a) shows a 
medium density system (with average distance between cross-links 
l c = 0.014//,), while the system in (b) shows a low density system 
(with l c = 0.019//,). Both have rigid cross-links, and average fila- 
ment length 0. 1 lip before removal of free end-points. 
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MC step/crosslink 

FIG. 7: Shear component h n , during a MC simulation under zero 
pressure, with and without slice moves. The simulated network is 
shown is Fig[6ta). As in Fig. [3] x and y denote directions of the two 
lattice vectors describing the geometry of this periodically replicated 
system. 



Our simulations focus explicitly on fluctuations in the posi- 
tions and orientations of cross-links, which primarily charac- 
terize the elasticity of this model system. In particular, a mi- 
crostate Y specifies only the configuration of cross-links, in- 
cluding the directions in which filaments pass through them. 
Thermal undulations of filament segments consistent with Y 
are integrated out beforehandaccording to the statistical me- 
chanics of a worm-like chain[18J. The "energy" E associ- 
ated with r thus in fact represents a free energy that accounts 
for the corresponding variety of chain configurations. It is 
a highly nonlinear function of cross-link arrangements, due 
to the inextensibility of a worm-like chain along its contour. 
These sharp nonlinearities foster heterogeneous stiffness and 
impede calculations of elastic response. 

Figs. |7] and |8] show results of Monte Carlo simulations for 
these model networks. They contrast fluctuations and relax- 
ation generated using conventional methods of sampling at 
constant pressure with those produced by slice moves. Slice 
shapes vy and Vy were chosen according to the recipe in 

Section HC] By setting v™ n = 0.05/zif l) , v™ in = 0.05^ t3 , 

v™ ax = 0.2h [ }f\ and v™ x = 0.2h^f\ where h""' is the shape 
matrix at the beginning of the simulation, we generate de- 
formation sub volumes with dimensions 5% — 20% of the ini- 
tial system size. The step size of trial shear deformations, 
e xy Ri 10 , was tuned during an equilibration phase of the 
simulation to establish an acceptance ratio of approximately 
0.5. 

Trajectories of spontaneous shear strain fluctuations are 
plotted in Fig [7] for the denser network configuration shown 
in Fig. |6ja). The enhanced efficiency offered by slice moves 
for sampling thermally accessible strain states is clearly evi- 
dent. By itself, the result for conventional, global strain moves 
provides no warning that it has failed to visit important re- 
gions of configuration space. One might therefore be tempted 
to estimate elastic susceptibilities, which would be orders of 
magnitude too small, from a severely deficient set of thermal 
fluctuations. 

The sparser network configuration shown in Fig[6j£>) is ex- 



0.25 r 




MC step/crosslink 

FIG. 8: Bulk component h xx during a MC simulation under applied 
pressure for a network undergoing collapse through buckling, with 
and without slice moves. The simulated network, whose horizontal 
direction is denoted by x, is shown in Fig[6j£>) . 



tremely susceptible to applied pressure. The bulk strain tra- 
jectory obtained using slice moves manifests this pliability 
through a systematic decrease in box size under load. (See 
Fig. [8]) With conventional methodology, by contrast, con- 
traction of the network as a whole necessitates deforming its 
densest regions, at least transiently. Indeed, when restricted 
to global strain moves, Monte Carlo sampling cannot access 
compressed states even within 10 7 sweeps. 

Because energy evaluations in network simulations are 
more numerically taxing than in hard disk simulations, the 
added overhead for performing slice moves amounts to a scant 
1% increase in run time compared to conventional simula- 
tions comprising the same total number of strain moves. This 
price is clearly outweighed by the dramatic gains in com- 
putational efficiency we have demonstrated. We expect effi- 
ciency considerations to similarly favor the use of slice moves 
for other complex systems that exhibit heterogeneous elastic- 
ity. In most physical contexts of interest, evaluating changes 
in potential energy due to intermolecular interactions will 
make negligible even the greatest expense brought on by slice 
moves, namely, determining which subvolume each particle 
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occupies when executing a deformation with non-rectangular 
slices. Furthermore, rectangular slices should suffice for ex- 
ploring many types of elasticity, e.g. in systems that are fluid; 
assigning particles to rectangular subvolumes is numerically 
inconsequential compared to calculating interaction energies 
for all but the simplest systems. 

IV. CONCLUSION 

We have shown how volume moves in constant-pressure 
simulations and strain moves in constant-stress simulations 
can be performed locally, such that intermolecular arrange- 
ments in much of a system remain undisturbed. Significant 
speedup of Monte Carlo simulations is expected for systems 
that are considerably nonuniform in stiffness. Example simu- 
lation results confirm that physically important strain states 
previously inaccessible as a matter of practice can now be 
readily explored. 

By facilitating spontaneous strain fluctuations, this method- 
ological advance promises to greatly extend the purview 
of techniques that assess linear elastic response via the 
fluctuation-dissipation theorem. Additionally, it provides a 
new type of collective Monte Carlo move as an alternative to 
cluster moves lll U Il9l 12011 . 

More broadly, it opens doors to applications in the many 
biophysical and materials contexts that involve spatially vary- 
ing density (as occurs in a material undergoing a phase transi- 
tion) and/or composition (as is routine in living cells). 
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APPENDIX A: ALGORITHM TO DETERMINE WHICH 
SLICE A POINT IS IN 



Executing a slice move requires determining the set of par- 
ticles that reside in each subvolume, before their coordinates 
can be appropriately transformed (according to Eq. [31 orlTsll. 
Performing this task efficiently is straightforward for subvol- 
umes that are rectangular in the reduced coordinate space. 
For non-rectangular slices, however, it can become both awk- 
ward and costly. Here we outline an algorithm that, for most 
points in a simulation box, reduces the classification problem 
to checking whether the point lies within a particular rectan- 
gle. 

The essence of this procedure is to inscribe a rectangle 
within each subvolume a (where a € {ii,v,vW,v^} in two 
dimensions). Particle coordinates can be quickly checked 
against these rectangles. Because useful deformation volumes 
tend to be small, most particles will fall within the inscribed 
rectangle of the undisturbed region u. Only a small fraction of 
particles need then be checked against subvolumes' full par- 
allelotope shapes. A systematic procedure for doing so is de- 
scribed below. 

Consider a particle located at position r in the reduced coor- 
dinate system, and a subvolume a centered at position c (also 
in the reduced coordinate system) with shape matrix oc !; -. We 
first determine which of the particle's periodic images, whose 
position we denote ?*, lies nearest c. We then compute a 
new set of reduced coordinates, r\ = (a -1 ),■_//! ; -^[r| — c#], refer- 
enced to the subvolume shape and translated so that the origin 
lies at c. If -1/2 < r\ < 1 /2 for all i = 1,2, ...,d, then the 
particle resides in a. By ordering subvolumes according to 
size, and checking particle positions against the largest slices 
first, we can ensure that most particles are assigned without 
numerous repetitions of these transformations. 
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